/************************************************************************/
/*  									*/
/*  lloyd.c								*/
/*									*/
/*  Tom Lookabaugh							*/
/*  August 13, 1987							*/
/*  Last Revision: August 14, 1987					*/
/*									*/
/*  DESCRIPTION								*/
/*	Runs the generalized Lloyd algorithm.  Attempts to replace 	*/
/*	empty cells by splitting the most populous cells.		*/
/*  RETURNS								*/
/*	Final distortion.  Modifies storage pointed to by cdbk, cdbkup	*/
/*	count, and trvect.						*/
/*  CALLS								*/
/*  	warning(), error(), distor(), cenupdt(), cenrplc(), perturb()	*/
/*									*/
/************************************************************************/

#include "stdvq.h"
#include <math.h>

double lloyd (cdbk, cdbkup, count, trvect, trvectin, cdbksz, threshold)

    double* cdbk;	/* codebook					*/
    double* cdbkup;	/* storage for updating codebook		*/
    long *count;	/* number of training sequence vectors per cell */
    DATATYPE *trvectin;	/* training sequence vector (input)		*/
    double *trvect;	/* training sequence vector			*/
    int cdbksz;		/* codebook size (in vectors)			*/
    double threshold;	/* convergence threshold			*/

{

    double dist;	/* distortion					*/
    int pass;		/* pass number					*/
    double olddist;	/* previous distortion				*/
    int emcltry;	/* attempt number for empty cell replacement	*/
    int c, d, j, k;	/* iteration counters				*/
    long t;		/* iteration counter				*/
    long emptcells;	/* number of empty cells			*/
    long *populcell;	/* indices of most populous cells		*/
    int diff, difftemp;	/* used in finding most populous cells		*/
    int replace;	/* cell to replace				*/
    int bestcdwd;	/* closest codeword to present data vector	*/
    

/* Loop: do while distortion is non-zero or fractional decrease in distortion
exceeds threshold.  Insure that at least two passes are completed.	*/

    dist = HUGE;
    pass = 1;
        
    do {
    
        olddist = dist;
        dist = 0.0;

/* Loop: repeat single Lloyd iteration upto EMCLTRIES number of times
(until there are no empty cells).					*/
    
	for (emcltry = 0; emcltry < EMCLTRIES; emcltry++) {
	
	    rewind(trsqfile);
	    for (c = 0; c < cdbksz; c++) {
	       for (d = 0; d < cdbkupdim; d++) cdbkup[c*cdbkupdim + d] = 0.0;
	       count[c] = 0;
	    }

/* Loop: read training sequence vectors, and perform full search through
codebook, updating as necessary for centroiding.			*/
	
	    for(t = 0; ; t++) {
	
	        if (fread(trvectin, sizeof(DATATYPE), dim, trsqfile) != dim)
	          break;
	        for (d = 0; d < dim; d++) trvect[d] = (double) trvectin[d];
	        dist += fullsearch(cdbk, cdbksz, trvect, &bestcdwd);
	        count[bestcdwd]++;
	        cenupdt(cdbkup + bestcdwd*cdbkupdim, trvect);
	    }
	
/* Centroid								*/


	    dist /= t;
	    emptcells = 0;
	    
	    for (c = 0; c < cdbksz; c++) {
                if (count[c] != 0) {
                    cenrplc(cdbkup + c*cdbkupdim, cdbk + c*cdbkdim, count[c]);
                }
                else emptcells++;
            }

/* If empty cells are present after codebook expansion,
         find highest population cells.                                 */
 
             if ((emptcells != 0) && (pass == 1)) {
         
                warning("empty cell replacement attempted (lloyd 1)");
                if ((populcell = (long *)calloc(emptcells, sizeof(long)))
                  == NULL) error("insufficient memory (lloyd 2)");
                k = 0;
                for (j = 0; ;j++) {
                    if (count[j] > 0) populcell[k++] = j;
                    if (k == emptcells) break;
                }
                if (k < emptcells) error("too many empty cells (lloyd 3)");
                for (j = emptcells; j < cdbksz; j++) {
                    diff = 0.0;
                    for (k = 0; k < emptcells; k++) {
                        if ((count[j] > 0) && ((difftemp = count[j]
                          - count[populcell[k]]) > diff)) {
                            diff = difftemp;
                            replace = k;
                        }
                    }
                    if (diff > 0) populcell[replace] = j;
                }

/* Replace empty cells with perturbed versions of high distortion cells.*/
            
                j = 0;
                for (k = 0; k < cdbksz; k++) {
                    if (count[k] == 0) {
                        perturb(cdbkup + cdbkupdim*populcell[j], cdbk + 
                          cdbkdim*k,
                          add_offset/(pow(2.0, (double) emcltry)),
                          mult_offset/(pow(2.0, (double) emcltry)));
                        j++;
                    }
                }
              
                free(populcell);
            }
         
            if (emptcells == 0) break;
        }
    
        if (emcltry == EMCLTRIES)
          error("couldn't eliminate empty cells (lloyd 4)");
        pass++;
        if (olddist < dist) warning("increase in distortion (lloyd 5)");
    }
    while ((dist != 0.0) && ((pass == 2) || ((olddist - dist)/dist
      > threshold) || (dist - olddist) > 0.0));
      
    return(dist);
}
 
	


